/// @file dag.cc
/// @brief Methods and helper classes in support of the Directed Acyclic Graph
///        data structure

// *************************************************************************

// Copyright© 2010 Texas A&M University, 
//                 College Station, Texas
//                 Contact: R Crosby <rwc@cs.tamu.edu>
//
// Copyright© 2003, 2004 by BiRC -- Bioinformatics Research Center
//                 University of Aarhus, Denmark
//                 Contact: Thomas Mailund <mailund@birc.dk>
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program.  If not, see <http://www.gnu.org/licenses/>.

// *************************************************************************

#include <algorithm>
#include <cassert>
#include <fstream>
#include <iomanip>
#include <string.h>

using namespace std;

#include "BuildFlags.hh"
#include "HashTable.hh"
#include "S.hh"
#include "UnivHashFixedVector.hh"
#include "dag.hh"
#include "options.hh"
#include "rootedtree.hh"

// *****************************************************************************
/// @brief A series of counters and methods in support of 
///        tracing the algorithm

namespace DagInfo {

    // Counters used during debugging and testing
    unsigned long long innerCountCalled       = 0;
    unsigned long long innerCountComputed     = 0;
    unsigned long long innerNodeCountCalled   = 0;
    unsigned long long innerNodeCountComputed = 0;
    unsigned long long targetTraversals       = 0;
    unsigned long long maxVectors             = 0;
    unsigned long long maxVectorMemory        = 0;
    unsigned long long maxVectorCapacity      = 0;
    unsigned long long nodeMemory             = 0; 
    unsigned long long innerNodes             = 0; 
    unsigned long long rootNodes              = 0; 
    unsigned long long leafNodes              = 0; 

    unsigned indent        = 0;		          // User to control trace indent

    const unsigned indentspace = 2;		  // Amount to increment indent

    inline void Entry ( const char *fn, void * obj ) {
	const char *lbl = "Entry: ";
	cerr << right << setw (strlen(lbl) + (indent++ * indentspace)) 
	     << lbl << fn << " @ " << obj << " ";
    }

    inline void Exit ( const char *fn, void * obj ) {
	const char *lbl = "Exit : ";
	cerr << right << setw (strlen(lbl) + (--indent * indentspace)) 
	     << lbl << fn << " @ " << obj << " ";
    }

    inline void DoFn ( const char *fn, void * obj ) {
	const char *lbl = "Do Fn: ";
	cerr << right << setw (strlen(lbl) + (indent * indentspace)) 
	     << lbl << fn << " @ " << obj << " ";
    }

}

#define ENFN(x) if (BuildFlags::trc) { DagInfo::Entry(x, this);
#define EXFN(x) if (BuildFlags::trc) { DagInfo::Exit (x, this);
#define DOFN(x) if (BuildFlags::trc) { DagInfo::DoFn (x, this);
#define ETRC() }
#define ENDL() endl << flush; ETRC()

// *****************************************************************************
/// @brief Base class for the various DAG entities

class DAGNode { 

    int          _a;            
    int          _b;
    int          _c;
    large_int_t  _nodecounted;
    CVEC         _dagcounted;
    int          _size;
    DAG &        _dag;

protected:

    /// @brief determine the color counts for a node
    void                SetColorCount(const S::colour_t c);

public:

    DAGNode     (DAG & d) : _a(0), _b(0), _c(0), _nodecounted(0), _size(0), _dag(d)
	{ }

    virtual             ~DAGNode     ()                                 {}

    //if nodecount has been called before this one will remember the outcome
    const large_int_t    nodecounted ()                              { return _nodecounted; }
    void                 nodecounted (const large_int_t nodecounted) { _nodecounted = nodecounted; }

    //if count has been called before this one will remember the outcome
    void                 dagcounted  (CVEC & dc)   { _dagcounted = dc; }
    CVEC &               dagcounted  ()            { return _dagcounted; }

    const int            size        ()            { return _size; }
    void                 size        (const int s) { _size = s; }

    /// Return a reference to the controlling DAG
    DAG &                Dag         ()            { return _dag; }
     
    const int  a ()             { return _a; }
    const int  b ()             { return _b; }
    const int  c ()             { return _c; }
    void       a (const int _a) { this->_a = _a; }
    void       b (const int _b) { this->_b = _b; }
    void       c (const int _c) { this->_c = _c; }

    void                 ResetVector    ();
    virtual void         ResetCount     () { _nodecounted = -1L; }

    virtual void         PruneVector    ( const unsigned nroots,
					  const unsigned rootIdx );

    virtual void         print          (ostream &os, const size_t indent) const;

    virtual void         colour_leaves  (S::colour_t c) = 0;

    /// @brief Return vector of counts

    /// Results[i] contains how many quartets is shared between the
    /// subtree rooted at this node and the tree given by root i of DAG
    virtual CVEC &       count          () = 0; 
		
    /// @brief Return single count value

    /// Returns how many quartets in the subtree is compatible with the
    /// current coloring
    virtual large_int_t  nodecount      () = 0;

    /// @brief Set the color counts for the node
    virtual void         SetColorCounts () = 0;

    /// @brief Used to obtain memory utilization information
    virtual unsigned     NodeMemory     () = 0;

    /// \brief Return the memory used by the active vectors
    virtual unsigned     VectorMemory   () { return sizeof (large_int_t) * dagcounted().size(); }

    /// \brief Return the memory allocated by the vectors
    virtual unsigned     VectorCapacity () { return sizeof (large_int_t) * dagcounted().capacity(); }
    
};

// *****************************************************************************
/// @brief A root node in the DAG

class DAGRoot : public DAGNode {

    DAGNode *      _child;
    S *            _s;
    const char *   _fname;			  // File name of original tree file
    const unsigned _lineno;			  // Line number within tree file

    /// @brief Number of taxa in this tree
    unsigned       _taxa;			  

public:
    DAGRoot       ( DAG &          d,
	            S *            s, 
		    const char*    fn,
		    const unsigned ln ) 
	: DAGNode(d), _s(s), _fname(fn), _lineno(ln), _taxa(1)
	{ }
    void        child         (DAGNode * c) 
	{ _child = c;
	    size(_child->size() + 1) ; }
    DAGNode*    child         () const { return _child; }
    void        print         (ostream &    os, 
			       const size_t indent) const;
    void        colour_leaves (S::colour_t  c);
    CVEC &      count         ();
    large_int_t nodecount     ();
    unsigned    NodeMemory    () { return sizeof ( DAGRoot ); }
    const char* Fname         () { return _fname; }
    unsigned    LineNo        () { return _lineno; }
    void        PruneVector   ( const unsigned nroots,
				const unsigned rootIdx) {}	// Don't prune the root

    void        SetColorCounts() { _child->SetColorCounts();
	                           SetColorCount(_s->Colour());
	                           a( _child->a() + a() );
	                           b( _child->b() + b() );
	                           c( _child->c() + c() );
                                  }

    /// @brief Accessor for the number of taxa in this tree
    unsigned    Taxa           () { return _taxa; }

    /// @brief Increment leaf count
    void        IncTaxa        () { _taxa++; }

};

// *****************************************************************************
/// @brief A inner node in the DAG

class DAGInner : public DAGNode {

    DAGNode *    _small;
    DAGNode *    _large;
    DAGRoot *    _lastRoot;

    /// \brief node count from the smaller subtree
    large_int_t  _s;
    /// \brief node count from the larger subtree
    large_int_t  _l;

    /// \brief Coefficients used in the dissimilar taxa case
    large_int_t  _a3_coef;
    large_int_t  _b3_coef;
    large_int_t  _c3_coef;
    large_int_t  _a32_coef;
    large_int_t  _b32_coef;
    large_int_t  _c32_coef;

public:
    DAGInner      ( DAG &     d,
		    DAGNode * small, 
		    DAGNode * large);

    DAGNode*     small         () const { return _small; }
    DAGNode*     large         () const { return _large; }

    void         print         ( ostream &    os, 
				 const size_t indent ) const;
    void         colour_leaves ( S::colour_t  c );
    CVEC &       count         ();
    large_int_t  nodecount     ();
    unsigned     NodeMemory    () { return sizeof ( DAGInner ); }

    void         LastRoot      ( DAGRoot * r ) { _lastRoot = r; }
    DAGRoot *    LastRoot      () { return _lastRoot; }

    void         SetColorCounts() { _small->SetColorCounts();
                                    _large->SetColorCounts();
 	                            a( _small->a() + _large->a() );
	                            b( _small->b() + _large->b() );
	                            c( _small->c() + _large->c() );
                                   }

    void         ResetCount     () { _s        = -1;
	                             _l        = -1;
				     _a3_coef  = -1;
				     _b3_coef  = -1;
				     _c3_coef  = -1;
				     _a32_coef = -1;
				     _b32_coef = -1;
				     _c32_coef = -1;
  	                             DAGNode::ResetCount(); }

};

// *****************************************************************************
/// @brief A leaf node in the DAG

class DAGLeaf : public DAGNode {

    S*  _s;

public:

    DAGLeaf       (DAG & d,
		   S *   s) 
	: DAGNode(d), _s(s) 
	{ size(1); }

    void        print         ( ostream &    s, 
			        const size_t indent ) const;
    void        colour_leaves ( S::colour_t  c );
    CVEC &      count         ();
    large_int_t nodecount     ();
    unsigned    NodeMemory    () { return sizeof ( DAGLeaf ); }

    void        SetColorCounts() { SetColorCount(_s->Colour()); }

};

// *****************************************************************************
/// @brief Hash table class for saving count vectors

class DAGHash 
    : public HashTable           < vector<S::colour_t>, 
				   CVEC >,
      public UnivHashFixedVector < S::colour_t >
{

    DAGHash( const unsigned nNodes,
	     const unsigned mapSize );

    virtual void     TableSize ( unsigned s );

    virtual unsigned KeyMem    ( const vector<S::colour_t> & v )
	{ return sizeof(S::colour_t) * v.size(); }

    virtual unsigned DataMem   ( const CVEC & v )
	{ return sizeof(large_int_t) * v.size(); }

    virtual unsigned DataCap   ( const CVEC & v )
	{ return sizeof(large_int_t) * v.capacity(); }

    friend class DAG;

};

// *****************************************************************************
/// @brief A visitor class to help in the merging of rooted tree entries
///        into the dag

class MergeIntoDAGVisitor : public RootedTreeVisitor {

    // Really the darkside of using visitors (or whatever to call them)
    // we have to handle the recursion our selves, thus this stack

    stack<DAGNode *> _dagnodes;
    DAG &            _dag;
    DAGRoot *        _currRoot;

public:
    MergeIntoDAGVisitor          (DAG  & dag) : _dag(dag) {}
    virtual ~MergeIntoDAGVisitor ()                       {}
    void     visit_root          (RootedTreeRoot *root,
				  const char *    fn,
				  const unsigned  ln );
    void     visit_leaf          (RootedTreeLeaf *leaf);
    void     visit_inner         (RootedTreeInner *inner);    

};

// *****************************************************************************
/// @brief Stream output routine to dump a vector of integers

ostream & operator<<(ostream &os, const vector<large_int_t> & results)
{ 

    os << "Vector[" << results.size() << "]: ";

    for ( unsigned i = 0; i < results.size(); i++ )
	os << results[i] << " ";
    return os; 
}

// *****************************************************************************
/// @brief Populate an inner node into the DAG

void
MergeIntoDAGVisitor::visit_inner(RootedTreeInner *inner)
{
    //recurse
    inner->left()->accept(this);
    inner->right()->accept(this);
    
    //update dag
    DAGNode *dn1 = _dagnodes.top(); _dagnodes.pop();
    DAGNode *dn2 = _dagnodes.top(); _dagnodes.pop();
    DAGInner *canonical = new DAGInner(_dag, dn1, dn2);
    pair<DAGNode*,DAGNode*> key = make_pair(canonical->small(),
					    canonical->large());

    map<pair<DAGNode *, DAGNode *>, DAGInner *>::iterator i = _dag._nodetable.find(key);
	
    if ( i != _dag._nodetable.end() ) {
	i->second->LastRoot(_currRoot);
	_dagnodes.push(_dag._nodetable[key]);
	delete canonical;
    }
    else {
	canonical->LastRoot(_currRoot);
	_dag._nodetable[key] = canonical;
	_dag._nodes.push_back(canonical);
	_dagnodes.push(canonical);
	DagInfo::innerNodes++;
    }    
}

// *****************************************************************************
/// @brief Populate a leaf node into the DAG

void
MergeIntoDAGVisitor::visit_leaf(RootedTreeLeaf *leaf)
{
    //update dag
    if (_dag._leaftable.count(leaf->s())) {
	_dagnodes.push(_dag._leaftable[leaf->s()]);
    }
    else {

	// If we found a new leaf after the first tree then set dissimilar mode.
	if ( _dag.Roots().size() > 0 )
	    _dag.Dissimilar(true);

	DAGLeaf *dl = new DAGLeaf(_dag, leaf->s());
	_dag._leaftable[leaf->s()] = dl;
	_dag._nodes.push_back(dl);
	_dagnodes.push(dl);
	DagInfo::leafNodes++;
    }  
    _currRoot->IncTaxa();
}

// *****************************************************************************
/// @brief Populate a root node into the DAG

void
MergeIntoDAGVisitor::visit_root(RootedTreeRoot *root,
				const char *    fn,
				const unsigned  ln )
{

    //FIXME if two identical trees are merged we should check it here
    //simply put a DAGNode in a table whenever we create a root
    _currRoot = new DAGRoot( _dag,
			     root->s(), 
			     fn, 
			     ln );

    //recurse
    root->child()->accept(this);
    _currRoot->child( _dagnodes.top() );

    _dagnodes.pop();
    _dag._roots.push_back(_currRoot);
    _dag._nodes.push_back(_currRoot);

    // If the number of taxa in this tree is not the same as the number of taxa
    // in the first tree set the disimilar flag.

    if ( _dag.Roots().size() > 1 &&		  // But not on the first root...
	 _currRoot->Taxa() != _dag.Roots()[0]->Taxa() )
	_dag.Dissimilar(true);

    DagInfo::rootNodes++;

}

// *****************************************************************************
/// @brief Basic destructor for the DAG
///
/// Prints the cache statistics if we're in verbose mode

DAG::~DAG()
{

    while ( !_nodes.empty() ) {
	delete _nodes.back();
	_nodes.pop_back();
    }

    if (options::verbose) {

	map<string, unsigned *> & stats = _hash->Stats();

	map<string, unsigned *>::const_iterator i;

	cerr << "Cache statistics:" << endl;
	for ( i = stats.begin(); i != stats.end(); i++ )
	    cerr << "  " 
		 << setw(14) << left << i->first << ": " 
		 << setw(11) << right << *(i->second) << endl;

    }
    

    delete _hash;

}

// ****************************************************************************
///
/// Loop through the set of roots building the count vectors for each one
/// at points along the way the vector set is pruned to the current length

void
DAG::count ( )
{

    PerfVar pf("DAGCount");
    PRF(pf, "Init");

    unsigned pruneCnt = Roots().size()/(100/_prunePct);

    unsigned rootIdx = 0;

    DagInfo::nodeMemory = NodeMemory();

    for ( _sourceIdx = 0; _sourceIdx < Roots().size(); _sourceIdx++) {

	ENFN("DAG::count") cerr << "Loop Start, Idx" << _sourceIdx << "> " << ENDL();

	if ( Dissimilar() )
	    _smap.ColourAll(S::None);

	DAGRoot * crt = Roots()[_sourceIdx];
	crt->colour_leaves(S::A);
	crt->count();

	EXFN("DAG::count") cerr << "Loop End, Idx<" << _sourceIdx << "> " << ENDL();
	DB() cerr << *this << endl;		  // Dump the dag at this point

	if (options::verbose) cerr << "." << flush;
	//if (options::verbose) PRF(pf, "Loop");

	if ( options::verbose && BuildFlags::perf ) 
	    // Get the vector usage information
	    VectorUsage ( DagInfo::maxVectorMemory,
			  DagInfo::maxVectorCapacity,
			  DagInfo::maxVectors );

	// Reset any vectors we're finished with
	DB() cerr << "Clearing <" << ClearList().size() << ">" << endl;
	while ( !ClearList().empty() ) {
	    ClearList().back()->ResetVector();
	    ClearList().pop_back();
	}

	// Time to make the donuts (prune the vectors)?
	if ( _sourceIdx && pruneCnt && ( _sourceIdx % pruneCnt == 0) ) {

	    if (options::verbose && BuildFlags::perf) cerr << endl;
	    PRF(pf, "Prune");
	    EXFN("DAG::count") cerr << "Prune" << ENDL();

	    for ( list<DAGNode*>::iterator i = Nodes().begin(); i != Nodes().end(); i++)

		if ( (*i)->dagcounted().size() > Roots().size() - rootIdx )

		    (*i)->PruneVector( Roots().size(),  rootIdx );
	    
	}

    }    

    if (options::verbose) cerr << endl;

}

// *****************************************************************************
/// @brief Return the file name associated with the specified index
/// 
/// @return Pointer to file name

const char *       
DAG::Fname ( const unsigned i ) 
{ 
    return _roots[i]->Fname(); 
}

// *****************************************************************************
/// @brief Return the line number associated with the specified index
/// 
/// @return Line number

const unsigned     
DAG::LineNo ( const unsigned i ) 
{  
    return _roots[i]->LineNo(); 
}

// *****************************************************************************
/// @brief Return the distribution of the vector lengths
/// 
/// @return Vector of lengths

vector<unsigned>
DAG::VSizeDist()
{

    vector<unsigned> v;
    v.resize(_roots.size());

    list<DAGNode *>::iterator i;
    for (i = _nodes.begin(); i != _nodes.end(); ++i) v[(*i)->dagcounted().size()]++;
    
    return v;

}

// *****************************************************************************
/// @brief Create the hash table
/// 
/// This was broken out of the constructor since the number of
/// nodes is required to allocate the table and that isn't 
/// known until the DAG is populated
///
/// @param mapSize Size of the SMap vector (key length)

void
DAG::MakeHashTable( const unsigned mapSize )
{

    _hash = new DAGHash ( _nodes.size(), mapSize );

}

// *****************************************************************************
/// @brief Clear the maps used to build the dag
/// 
/// This can't be done until all the trees are loaded

void
DAG::ClearMaps()
{

    _nodetable.clear();
    _leaftable.clear();

}


// *****************************************************************************
/// @brief Merge a rooted tree into the dag
/// 
/// Basically a postorder traversal of the tree
///
/// @param rt Reference to a rooted tree to merge into the dag
/// @param fn Tree file name
/// @param rt Line number of tree in input file

void
DAG::merge ( RootedTree &   rt, 
	     const char *   fn,
	     const unsigned ln )
{

    MergeIntoDAGVisitor dv(*this);		  // This guarantees that it's destroyed
    rt.root()->accept( &dv, fn, ln );

    _labelWidth = max( _labelWidth, (unsigned)strlen(fn) );
    
}

// *****************************************************************************
/// Counts how many quartets are compatible with the current
/// coloring. Results[i] contains the number of quartets in tree i
/// compatible with the current coloring
/// 
/// @param results Reference to the vector of counts to populate

void
DAG::nodecount(CVEC & results)
{

    ENFN("DAG::nodecount") cerr << results << ENDL();

    ResetCount();

    // Loop through all subsequent roots

    for ( _targetIdx = _sourceIdx + 1; _targetIdx < _roots.size(); _targetIdx++) {

	if ( Dissimilar() )
	    _roots[_targetIdx]->SetColorCounts();

	results.push_back(_roots[_targetIdx]->nodecount());
    }

    EXFN("DAG::nodecount") cerr << results << ENDL();

}

// *****************************************************************************
/// @brief Return DAG memory usage
///
/// @return Memory (bytes) used by dag
///

unsigned long long
DAG::NodeMemory ()
{
    unsigned long long mz = 0;

    for ( list<DAGNode*>::const_iterator i = Nodes().begin(); i != Nodes().end(); i++ )
	mz += (*i)->NodeMemory();

    return mz;
}

// *****************************************************************************
/// @brief Return DAG vector usage information
///
/// @param maxAlloc Reference to high water mark memory for used vectors
/// @param maxCap   Reference to high water mark for pre-allocated capacity
/// @param maxVec   Reference to number of active vectors at maxAlloc point

void 
DAG::VectorUsage ( unsigned long long & maxAlloc,
		   unsigned long long & maxCap,
		   unsigned long long & maxVec )
{
    unsigned long long ma = 0;
    unsigned long long mc = 0;
    unsigned           vc = 0;

    for ( list<DAGNode*>::const_iterator i = Nodes().begin(); i != Nodes().end(); i++ ) {

	if ( ! (*i)->dagcounted().empty() ) {
	    vc++;
	    ma += (*i)->VectorMemory();
	}

	mc += (*i)->VectorCapacity();
    }

    if ( ma > maxAlloc ) {
	maxAlloc = ma;
	maxVec   = vc;
    }

    if ( mc > maxCap ) 
	maxCap = mc;

}

// *****************************************************************************
/// @brief Stream helper function to output the DAG
///
/// @param os Output stream reference

void
DAG::print(ostream &os) const
{
 
    {
	vector<DAGRoot *>::const_iterator i;
	os << "DAG roots {" << endl;
	for (i = _roots.begin(); i != _roots.end(); ++i) (*i)->print(os, 2);
	os << "}" << endl;
    }

    // {
    //  	list<DAGNode *>::const_iterator i;
    //  	os << "DAG nodes {" << endl;
    //  	for (i = _nodes.begin(); i != _nodes.end(); ++i) (*i)->print(os, 2);
    //  	os << "}" << endl;
    //  }

}

// *****************************************************************************
/// @brief Zero the nodecount for the entire DAG
///
/// Iterates through the entire DAG zeroing the nodecount values.
/// 

void
DAG::ResetCount()
{

    DOFN("DAG::ResetCount") cerr << ENDL();

    for( list<DAGNode *>::iterator i = _nodes.begin(); i != _nodes.end(); ++i)
	(*i)->ResetCount();

}

// *****************************************************************************
/// @brief Clear the dagcounted vector for the entire DAG
///
/// This iterates through all the nodes in the DAG clearing the dagcounted 
/// vector 
/// 

void
DAG::ResetVector()
{

    DOFN("DAG::ResetVector") cerr << ENDL();

    for( list<DAGNode *>::iterator i = _nodes.begin(); i != _nodes.end(); ++i)
	(*i)->ResetVector();

}

// *****************************************************************************
/// Return the number of taza in the specified tree
/// Not inlined since DAGRoot is not available outside of this source file

unsigned
DAG::Taxa (const unsigned tree)
{
    return Roots()[tree]->Taxa();
}
// *****************************************************************************
/// @brief Return the count vector for the specified node
///
/// This is not inlined into the .hh because it needs access to the dagcount()
/// method of DAGNode.

CVEC &
DAG::operator[] ( unsigned i )
{

    return Roots()[i]->dagcounted();

}

// *****************************************************************************
/// @brief Create a hash table 
///
/// Use the number of DAG nodes to guess the number of slots in the hash table.
/// The size of the S color vector is passed as the hash key vector size.
///
/// @param nNodes Number of DAG nodes 
/// @param mapSize Size of the S vector of colours

DAGHash::DAGHash( const unsigned nNodes,
		  const unsigned mapSize )
    : UnivHashFixedVector< S::colour_t >(mapSize)
{

    // Tell hash table about the encoder - this needs to happen first!
    Encoder(dynamic_cast < UnivHashFixedVector < S::colour_t > * >( this ));

    // Set hash table size
    TableSize(nNodes);

}

// *****************************************************************************
/// @brief Set the size of the hash table
///
/// Use the number of DAG nodes to guess the table size.
/// The value passed to the hash table constructor is used to
/// call the virtual table size method which ends up here...
///
/// @param s Value passed to constructor

void
DAGHash::TableSize ( unsigned s ) 
{

    /// @todo Come up with a reasonable algorithm for setting table size

    s = 14;

    DB() cerr << "DAGHash::TableSize s<" << s << ">" << endl << flush;

    HashTable< vector<S::colour_t>, vector<large_int_t> >::TableSize(s); 


}


// *****************************************************************************
/// @brief Constructor for the inner nodes
///
/// Construct an inner node with a large and small subtree
///
/// @param small Pointer to the smaller subtree
/// @param large Pointer to the larger subtree

DAGInner::DAGInner(DAG & d, DAGNode *small, DAGNode *large) :
    DAGNode(d), 
    _s(-1), _l(-1),
    _a3_coef(-1),  _b3_coef(-1),  _c3_coef(-1),
    _a32_coef(-1), _b32_coef(-1), _c32_coef(-1)
{
 
    if(small->size() > large->size()) {
	_small = large;
	_large = small;
    } 
    else {
	_small = small;
	_large = large;
    }

    size(_small->size() + _large->size() + 1);

}

// *****************************************************************************
/// @brief Color all leaves below current inner node
///
/// @param c Color to set

void
DAGInner::colour_leaves(S::colour_t c)
{

    ENFN("DAGInner::colour_leaves") cerr << c << ENDL();

    small()->colour_leaves(c);
    large()->colour_leaves(c);

    EXFN("DAGLeaf::colour_leaves") cerr << c << ENDL();

}

// *****************************************************************************
/// @brief Generate count vector for an inner node
///
/// @return Vector of counts

CVEC &
DAGInner::count()
{

    ENFN("DAGInner::count") cerr << ENDL();

    // If we're done with this inner node, mark it to have it's vector cleared
    if ( _lastRoot == Dag().Roots()[Dag().SourceIdx()]) 
	Dag().ClearList().push_back(this);

    DagInfo::innerCountCalled++;

    if( !dagcounted().empty() ) {

	colour_leaves(S::C);
	EXFN("DAGInner::count") cerr << "Already Determined: " << dagcounted() << ENDL();	
	return dagcounted();
    }
  
    DagInfo::innerCountComputed++;

    small()->colour_leaves(S::B);
    DOFN("DAGInner::count") cerr << "set small(B) " ; Dag().Smap().PrintTable(cerr); cerr << ENDL();

    // CVEC * xp = Dag().Dissimilar() ? NULL
    //                                : Dag().Hash()->Find ( Dag().Smap().ColorVector() );
    CVEC * xp = Dag().Hash()->Find ( Dag().Smap().ColorVector() );
    if (!xp ) {
	xp = new CVEC;
	Dag().nodecount(*xp);
	Dag().Hash()->Insert( Dag().Smap().ColorVector(), *xp );
	DOFN("DAGInner::count") cerr << "NotFound " << ENDL();
    }
    else {
	DOFN("DAGInner::count") cerr << "Found " << ENDL();
    }

    small()->colour_leaves(S::C);
    DOFN("DAGInner::count") cerr << "set small(C) " ; Dag().Smap().PrintTable(cerr); cerr << ENDL();
    CVEC & y = large()->count(); 

    small()->colour_leaves(S::A);
    DOFN("DAGInner::count") cerr << "set small(A) " ; Dag().Smap().PrintTable(cerr); cerr << ENDL();
    CVEC & z = small()->count();

    // The array may be of different sizes due to the timing of their
    // creation. The offsets are used to make sure things align correctly
    unsigned int iMin = min ( xp->size(), 
			      min ( y.size(), 
				    min ( z.size(), Dag().Roots().size() - Dag().SourceIdx() - 1 ))) ;

    unsigned int xOff = xp->size() - iMin;
    unsigned int yOff = y.size() - iMin;
    unsigned int zOff = z.size() - iMin;

    // The size of the array is based on the number of slots to the right of the
    // current pivot as representede by the current root index
    dagcounted().reserve(Dag().Roots().size() - Dag().SourceIdx() - 1);

    for(unsigned int i = 0; i < iMin; ++i)
	dagcounted().push_back ( xp->at(i + xOff) + y.at(i + yOff) + z.at(i + zOff) );

    delete xp;

    EXFN("DAGInner::count") cerr << "Computed:  " << dagcounted();
    Dag().Smap().PrintTable(cerr); 
    cerr << ENDL();	

    return dagcounted();    
}

// *****************************************************************************
/// @brief Determine the count for an inner node
///
/// This routine contains the guts of the counting process computing
/// the count based on the colors of the smaller and larger subtrees
///
/// @return Number of matching nodes
 
large_int_t
DAGInner::nodecount()
{

    ENFN("DAGInner::nodecount") cerr << ENDL();

    DagInfo::innerNodeCountCalled++;

    //if this one has been nodecounted before return that
    //a speedup in the naive approach, won't help us in the bottom up
    //construction if that is worth doing

    if ( !Dag().Dissimilar() && nodecounted() != -1L ) {
	EXFN("DAGInner::nodecount") cerr << "Already Determined " << nodecounted() << ENDL();
	return nodecounted();
    }
    
    // The following needs to happen before the rest of the computations in order to set the color counts
    // in the similar taxa case
    if ( _s == -1 ) {
	_s = small()->nodecount();
	_l = large()->nodecount();
	DagInfo::innerNodeCountComputed++;
    } 

    DOFN("DAGInner::nodecount") cerr << "Color Counts (" << a() << ',' << b() << ',' << c() << ")" << ENDL();

    large_int_t q;

    if ( !Dag().Dissimilar() ) {

	a( small()->a() + large()->a() );
	b( small()->b() + large()->b() );
	c( small()->c() + large()->c() );

	//these used to be ints
	large_int_t a1 = (unsigned long)small()->a();
	large_int_t b1 = (unsigned long)small()->b();
	large_int_t c1 = (unsigned long)small()->c();
	large_int_t a2 = (unsigned long)large()->a();
	large_int_t b2 = (unsigned long)large()->b();
	large_int_t c2 = (unsigned long)large()->c();

	large_int_t a3 = (unsigned long)Dag().Smap().Numberof(S::A) - a1 - a2;
	large_int_t b3 = (unsigned long)Dag().Smap().Numberof(S::B) - b1 - b2;
	large_int_t c3 = (unsigned long)Dag().Smap().Numberof(S::C) - c1 - c2;

	DOFN("DAGInner::nodecount") cerr << "S/M/L values (" << a1 << ',' << b1 << ',' << c1 << ") (" << a2<< ','  << b2 << ',' << c2 << ") (" << a3 << ',' << b3 << ',' << c3 << ")" << ENDL();

	//
	// The two sets of equations are equivalent
	// the first is used to provide a slightly better performance profile
	// 

	q = ( 
	    (a3*a3) * ( b2 * c1 + b1 * c2 ) + 
	    a3 * ( b2 * c1 * (b2 + c1 - 3) + b1 * c2 * (b1 + c2 - 3) ) +
	    (a2*a2) * (b3 * c1 + b1 * c3) + 
	    a1 * (b3 * c2 * ( a1 + b3 + c2 - 3) + b2 * (a1 + b2 - 3) * c3 + b2 * (c3*c3)) +
	    a2 * (b3 * c1 * (b3 + c1 - 3) + b1 * c3 * (b1 + c3 - 3))
	    ) / 2;

    }
    else {

	//these used to be ints
	large_int_t a1 = (unsigned long)small()->a();
	large_int_t b1 = (unsigned long)small()->b();
	large_int_t c1 = (unsigned long)small()->c();
	large_int_t a2 = (unsigned long)large()->a();
	large_int_t b2 = (unsigned long)large()->b();
	large_int_t c2 = (unsigned long)large()->c();

	large_int_t a3 = (unsigned long)Dag().Roots()[Dag().TargetIdx()]->a() - a1 - a2;
	large_int_t b3 = (unsigned long)Dag().Roots()[Dag().TargetIdx()]->b() - b1 - b2;
	large_int_t c3 = (unsigned long)Dag().Roots()[Dag().TargetIdx()]->c() - c1 - c2;

	DOFN("DAGInner::nodecount") cerr << "S/M/L values (" << a1 << ',' << b1 << ',' << c1 << ") (" << a2<< ','  << b2 << ',' << c2 << ") (" << a3 << ',' << b3 << ',' << c3 << ")" << ENDL();

	if ( _a3_coef == -1 ) {
	    _a3_coef  = b2 * c1 * (b2 + c1 - 3) + b1 * c2 * (b1 + c2 - 3);
	    _a32_coef = b2 * c1                 + b1 * c2;
	    _b3_coef  = a2 * c1 * (a2 + c1 - 3) + a1 * c2 * (a1 + c2 - 3);
	    _b32_coef = a2 * c1                 + a1 * c2;
	    _c3_coef  = a2 * b1 * (a2 + b1 - 3) + a1 * b2 * (a1 + b2 - 3);
	    _c32_coef = a2 * b1                 + a1 * b2;
	}

	//
	// The two sets of equations are equivalent
	// the first is used to provide a slightly better performance profile
	// 

	q = ( 
	    a3 * _a3_coef + a3 * a3 * _a32_coef +
	    b3 * _b3_coef + b3 * b3 * _b32_coef +
	    c3 * _c3_coef + c3 * c3 * _c32_coef
	    ) / 2;

    }

    // large_int_t q = 0UL;
    // q = q + (a1*(a1-1)*(b2*c3+b3*c2))/2;
    // q = q + (a2*(a2-1)*(b1*c3+b3*c1))/2;
    // q = q + (a3*(a3-1)*(b1*c2+b2*c1))/2;
    // q = q + (b1*(b1-1)*(a2*c3+a3*c2))/2;
    // q = q + (b2*(b2-1)*(a1*c3+a3*c1))/2;
    // q = q + (b3*(b3-1)*(a1*c2+a2*c1))/2;
    // q = q + (c1*(c1-1)*(a2*b3+a3*b2))/2;
    // q = q + (c2*(c2-1)*(a1*b3+a3*b1))/2;
    // q = q + (c3*(c3-1)*(a1*b2+a2*b1))/2;

    nodecounted(_s + _l + q);

    EXFN("DAGInner::nodecount") cerr << "Computed " << nodecounted() 
				     << " A/B/C<" << Dag().Smap().Numberof(S::A) 
				     << "," << Dag().Smap().Numberof(S::B)
				     << "," << Dag().Smap().Numberof(S::C) << ">" << ENDL();

    return nodecounted();
}

// *****************************************************************************
/// @brief Stream helper method to print the contents of an inner node
///
/// @param os Reference to the output stream to use
/// @param indent Amount to indent the printed output

void
DAGInner::print(ostream &os, const size_t indent) const
{
    {
	const char *lbl = "Inner: ";
	os << right << setw (strlen(lbl) + indent) << lbl << this;
	DAGNode::print(os, indent);
	os << endl;
    }

    {
	const char *lbl = "-> small: ";
//	os << right << setw (strlen(lbl) + indent) << lbl << small() << endl;
	os << right << setw (strlen(lbl) + indent) << lbl;
	small()->print(os, indent + 2);
    }

    {
	const char *lbl = "-> large: ";
//	os << right << setw (strlen(lbl) + indent)<< lbl << large() << endl;
	os << right << setw (strlen(lbl) + indent)<< lbl;
	large()->print(os, indent + 2);
    }


}

// *****************************************************************************
/// @brief Update the color of the node
///
/// @param c Color to set on the leaf

void
DAGLeaf::colour_leaves(S::colour_t c)
{

    DOFN("DAGLeaf::colour_leaves") cerr << *(this->_s) << " New Color: " << c <<ENDL();

    _s->Colour(c);

}

// *****************************************************************************
/// @brief Set the color to C and return the count vector
///
/// @return Vector of count pointers

CVEC &
DAGLeaf::count()
{

    ENFN("DAGLeaf::count") cerr << *(this->_s) << ENDL();

    colour_leaves(S::C);

    // 
    // The size of the array is based on the number of slots to the right of the
    // current pivot as representede by the current root index
    // 

    dagcounted().resize(Dag().Roots().size() - Dag().SourceIdx() - 1, 0);

    EXFN("DAGLeaf::count") cerr << dagcounted() << ENDL();	

    return dagcounted();

}

// *****************************************************************************
/// @brief Set the initial count for the node based on color
///
/// @return Count (zero)

large_int_t
DAGLeaf::nodecount()
{

    DOFN("DAGLeaf::nodecount") cerr << *(this->_s) << ENDL();

    if ( !Dag().Dissimilar() )
	SetColorCount(_s->Colour());

    return 0;

}

// *****************************************************************************
/// @brief Stream helper function to print the leaf node
///
/// @param os Reference to the output stream to use
/// @param indent Quantity to indent the output

void
DAGLeaf::print( ostream &os, 
	        const size_t indent ) const
{
    const char *lbl = "Leaf:  ";
    os << right << setw (strlen(lbl) + indent) << lbl << this << " " << *_s;
    DAGNode::print(os, indent);
    os << endl;

}

// *****************************************************************************
/// @brief Helper function to output the contents of a DAG node
///
/// @param os Reference to an output stream to write to
/// @param indent amount to indent the printed output

void
DAGNode::print ( ostream& os, 
		 const size_t  indent) const
{

    os << " a/b/c<" << _a << "," << _b << "," << _c << ">";
    os << " nodecounted<" << _nodecounted << ">";
    os << " size<" << _size << ">";
    os << " dagcounted " << _dagcounted ;

}

// *****************************************************************************
/// @brief Trim the first rows off based on the current root index
///
/// @param nroots Total number of rows
/// @param rootIdx Current index into root table

void
DAGNode::PruneVector ( const unsigned nroots,
		       const unsigned rootIdx )
{

    CVEC nv;

    nv.resize(nroots - rootIdx);

    copy ( _dagcounted.begin() + rootIdx - nroots + _dagcounted.size(), 
	   _dagcounted.end(), 
	   nv.begin() );

    _dagcounted.swap(nv);

}

// *****************************************************************************
/// @brief Reset the count vector

void
DAGNode::ResetVector() 
{
    CVEC v;
    dagcounted().swap(v);
}

// *****************************************************************************
/// @brief Set the color counters for the leaf

void
DAGNode::SetColorCount(const S::colour_t co)
{

    switch (co)  {

    case S::A:
	a(1); b(0); c(0); 
	break;
    case S::B:
	a(0); b(1); c(0); 
	break;
    case S::C:
	a(0); b(0); c(1); 
	break;
    case S::None:
	a(0); b(0); c(0); 
	break;
    default:
	assert(0);

    }

    DOFN("DAGNode::SetColorCount") cerr << "Color Counts (" << a() << ',' << b() << ',' << c() << ")" << ENDL();

}

// *****************************************************************************
/// @brief This sets the color of the root node
///
/// @param c Color

void
DAGRoot::colour_leaves(S::colour_t c)
{

    DOFN("DAGRoot::colour_leaves") cerr << *(this->_s) << " New Color: " << c <<ENDL();
    
    _s->Colour(c);
    _child->colour_leaves(c);

}

// *****************************************************************************
/// @brief Build a vector of counts for a root element
///
/// Color the root only according to "C" and then count the children
///
/// @return Vector of counts

CVEC &
DAGRoot::count()
{

    ENFN("DAGRoot::count") cerr << ENDL();

    _s->Colour(S::C);

    dagcounted() = child()->count();

    EXFN("DAGRoot::count") cerr << ENDL();
    return dagcounted();

}

// *****************************************************************************
/// @brief Get the count for a tree
///
/// Count the children
///
/// @return Count

large_int_t
DAGRoot::nodecount()
{

    ENFN("DAGRoot::nodecount") cerr  << ENDL();
    
    DagInfo::targetTraversals++;

    large_int_t rv = child()->nodecount();

    EXFN("DAGRoot::nodecount") cerr  << ENDL();

    return rv;
}

// *****************************************************************************
/// @brief Print the contents of the dag root element
///
/// Used as a helper function to the stream output for DAGRoot
///
/// @param os Output stream
/// @param indent Amount to indent printed output

void
DAGRoot::print(ostream &os, const size_t indent) const
{
    const char *lbl = "Root:  ";
    os << right << setw (strlen(lbl) + indent) << lbl << this << " " << *_s;

    DAGNode::print(os, indent);

    os << " <" << _fname << "," << _lineno << ">";

    os << " Taxa<" << _taxa << ">";

    os << endl;

    child()->print(os, indent + 2);

}
 
